home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / complex / math.c next >
Encoding:
C/C++ Source or Header  |  2001-11-01  |  21.1 KB  |  1,000 lines
  1. /* complex/math.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jorma Olavi TΣhtinen, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Basic complex arithmetic functions
  21.  
  22.  * Original version by Jorma Olavi TΣhtinen <jotahtin@cc.hut.fi>
  23.  *
  24.  * Modified for GSL by Brian Gough, 3/2000
  25.  */
  26.  
  27. /* The following references describe the methods used in these
  28.  * functions,
  29.  *
  30.  *   T. E. Hull and Thomas F. Fairgrieve and Ping Tak Peter Tang,
  31.  *   "Implementing Complex Elementary Functions Using Exception
  32.  *   Handling", ACM Transactions on Mathematical Software, Volume 20
  33.  *   (1994), pp 215-244, Corrigenda, p553
  34.  *
  35.  *   Hull et al, "Implementing the complex arcsin and arccosine
  36.  *   functions using exception handling", ACM Transactions on
  37.  *   Mathematical Software, Volume 23 (1997) pp 299-335
  38.  *
  39.  *   Abramowitz and Stegun, Handbook of Mathematical Functions, "Inverse
  40.  *   Circular Functions in Terms of Real and Imaginary Parts", Formulas
  41.  *   4.4.37, 4.4.38, 4.4.39
  42.  */
  43.  
  44. #include <config.h>
  45. #include <math.h>
  46. #include <gsl/gsl_math.h>
  47. #include <gsl/gsl_complex.h>
  48. #include <gsl/gsl_complex_math.h>
  49.  
  50. /**********************************************************************
  51.  * Complex numbers
  52.  **********************************************************************/
  53.  
  54. #ifndef HIDE_INLINE_STATIC
  55. gsl_complex
  56. gsl_complex_rect (double x, double y)
  57. {                /* return z = x + i y */
  58.   gsl_complex z;
  59.   GSL_SET_COMPLEX (&z, x, y);
  60.   return z;
  61. }
  62. #endif
  63.  
  64. gsl_complex
  65. gsl_complex_polar (double r, double theta)
  66. {                /* return z = r exp(i theta) */
  67.   gsl_complex z;
  68.   GSL_SET_COMPLEX (&z, r * cos (theta), r * sin (theta));
  69.   return z;
  70. }
  71.  
  72. /**********************************************************************
  73.  * Properties of complex numbers
  74.  **********************************************************************/
  75.  
  76. double
  77. gsl_complex_arg (gsl_complex z)
  78. {                /* return arg(z),  -pi < arg(z) <= +pi */
  79.   return atan2 (GSL_IMAG (z), GSL_REAL (z));
  80. }
  81.  
  82. double
  83. gsl_complex_abs (gsl_complex z)
  84. {                /* return |z| */
  85.   return hypot (GSL_REAL (z), GSL_IMAG (z));
  86. }
  87.  
  88. double
  89. gsl_complex_abs2 (gsl_complex z)
  90. {                /* return |z|^2 */
  91.   double x = GSL_REAL (z);
  92.   double y = GSL_IMAG (z);
  93.  
  94.   return (x * x + y * y);
  95. }
  96.  
  97. double
  98. gsl_complex_logabs (gsl_complex z)
  99. {                /* return log|z| */
  100.   double xabs = fabs (GSL_REAL (z));
  101.   double yabs = fabs (GSL_IMAG (z));
  102.   double max, u;
  103.  
  104.   if (xabs >= yabs)
  105.     {
  106.       max = xabs;
  107.       u = yabs / xabs;
  108.     }
  109.   else
  110.     {
  111.       max = yabs;
  112.       u = xabs / yabs;
  113.     }
  114.  
  115.   /* Handle underflow when u is close to 0 */
  116.  
  117.   return log (max) + 0.5 * log1p (u * u);
  118. }
  119.  
  120.  
  121. /***********************************************************************
  122.  * Complex arithmetic operators
  123.  ***********************************************************************/
  124.  
  125. gsl_complex
  126. gsl_complex_add (gsl_complex a, gsl_complex b)
  127. {                /* z=a+b */
  128.   double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  129.   double br = GSL_REAL (b), bi = GSL_IMAG (b);
  130.  
  131.   gsl_complex z;
  132.   GSL_SET_COMPLEX (&z, ar + br, ai + bi);
  133.   return z;
  134. }
  135.  
  136. gsl_complex
  137. gsl_complex_add_real (gsl_complex a, double x)
  138. {                /* z=a+x */
  139.   gsl_complex z;
  140.   GSL_SET_COMPLEX (&z, GSL_REAL (a) + x, GSL_IMAG (a));
  141.   return z;
  142. }
  143.  
  144. gsl_complex
  145. gsl_complex_add_imag (gsl_complex a, double y)
  146. {                /* z=a+iy */
  147.   gsl_complex z;
  148.   GSL_SET_COMPLEX (&z, GSL_REAL (a), GSL_IMAG (a) + y);
  149.   return z;
  150. }
  151.  
  152.  
  153. gsl_complex
  154. gsl_complex_sub (gsl_complex a, gsl_complex b)
  155. {                /* z=a-b */
  156.   double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  157.   double br = GSL_REAL (b), bi = GSL_IMAG (b);
  158.  
  159.   gsl_complex z;
  160.   GSL_SET_COMPLEX (&z, ar - br, ai - bi);
  161.   return z;
  162. }
  163.  
  164. gsl_complex
  165. gsl_complex_sub_real (gsl_complex a, double x)
  166. {                /* z=a-x */
  167.   gsl_complex z;
  168.   GSL_SET_COMPLEX (&z, GSL_REAL (a) - x, GSL_IMAG (a));
  169.   return z;
  170. }
  171.  
  172. gsl_complex
  173. gsl_complex_sub_imag (gsl_complex a, double y)
  174. {                /* z=a-iy */
  175.   gsl_complex z;
  176.   GSL_SET_COMPLEX (&z, GSL_REAL (a), GSL_IMAG (a) - y);
  177.   return z;
  178. }
  179.  
  180. gsl_complex
  181. gsl_complex_mul (gsl_complex a, gsl_complex b)
  182. {                /* z=a*b */
  183.   double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  184.   double br = GSL_REAL (b), bi = GSL_IMAG (b);
  185.  
  186.   gsl_complex z;
  187.   GSL_SET_COMPLEX (&z, ar * br - ai * bi, ar * bi + ai * br);
  188.   return z;
  189. }
  190.  
  191. gsl_complex
  192. gsl_complex_mul_real (gsl_complex a, double x)
  193. {                /* z=a*x */
  194.   gsl_complex z;
  195.   GSL_SET_COMPLEX (&z, x * GSL_REAL (a), x * GSL_IMAG (a));
  196.   return z;
  197. }
  198.  
  199. gsl_complex
  200. gsl_complex_mul_imag (gsl_complex a, double y)
  201. {                /* z=a*iy */
  202.   gsl_complex z;
  203.   GSL_SET_COMPLEX (&z, -y * GSL_IMAG (a), y * GSL_REAL (a));
  204.   return z;
  205. }
  206.  
  207. gsl_complex
  208. gsl_complex_div (gsl_complex a, gsl_complex b)
  209. {                /* z=a/b */
  210.   double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  211.   double br = GSL_REAL (b), bi = GSL_IMAG (b);
  212.  
  213.   double s = 1.0 / gsl_complex_abs (b);
  214.  
  215.   double sbr = s * br;
  216.   double sbi = s * bi;
  217.  
  218.   double zr = (ar * sbr + ai * sbi) * s;
  219.   double zi = (ai * sbr - ar * sbi) * s;
  220.  
  221.   gsl_complex z;
  222.   GSL_SET_COMPLEX (&z, zr, zi);
  223.   return z;
  224. }
  225.  
  226. gsl_complex
  227. gsl_complex_div_real (gsl_complex a, double x)
  228. {                /* z=a/x */
  229.   gsl_complex z;
  230.   GSL_SET_COMPLEX (&z, GSL_REAL (a) / x, GSL_IMAG (a) / x);
  231.   return z;
  232. }
  233.  
  234. gsl_complex
  235. gsl_complex_div_imag (gsl_complex a, double y)
  236. {                /* z=a/(iy) */
  237.   gsl_complex z;
  238.   GSL_SET_COMPLEX (&z, GSL_IMAG (a) / y,  - GSL_REAL (a) / y);
  239.   return z;
  240. }
  241.  
  242. gsl_complex
  243. gsl_complex_conjugate (gsl_complex a)
  244. {                /* z=conj(a) */
  245.   gsl_complex z;
  246.   GSL_SET_COMPLEX (&z, GSL_REAL (a), -GSL_IMAG (a));
  247.   return z;
  248. }
  249.  
  250. gsl_complex
  251. gsl_complex_negative (gsl_complex a)
  252. {                /* z=-a */
  253.   gsl_complex z;
  254.   GSL_SET_COMPLEX (&z, -GSL_REAL (a), -GSL_IMAG (a));
  255.   return z;
  256. }
  257.  
  258. gsl_complex
  259. gsl_complex_inverse (gsl_complex a)
  260. {                /* z=1/a */
  261.   double s = 1.0 / gsl_complex_abs (a);
  262.  
  263.   gsl_complex z;
  264.   GSL_SET_COMPLEX (&z, (GSL_REAL (a) * s) * s, -(GSL_IMAG (a) * s) * s);
  265.   return z;
  266. }
  267.  
  268. /**********************************************************************
  269.  * Elementary complex functions
  270.  **********************************************************************/
  271.  
  272. gsl_complex
  273. gsl_complex_sqrt (gsl_complex a)
  274. {                /* z=sqrt(a) */
  275.   gsl_complex z;
  276.  
  277.   if (GSL_REAL (a) == 0.0 && GSL_IMAG (a) == 0.0)
  278.     {
  279.       GSL_SET_COMPLEX (&z, 0, 0);
  280.     }
  281.   else
  282.     {
  283.       double x = fabs (GSL_REAL (a));
  284.       double y = fabs (GSL_IMAG (a));
  285.       double w;
  286.  
  287.       if (x >= y)
  288.     {
  289.       double t = y / x;
  290.       w = sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 + t * t)));
  291.     }
  292.       else
  293.     {
  294.       double t = x / y;
  295.       w = sqrt (y) * sqrt (0.5 * (t + sqrt (1.0 + t * t)));
  296.     }
  297.  
  298.       if (GSL_REAL (a) >= 0.0)
  299.     {
  300.       double ai = GSL_IMAG (a);
  301.       GSL_SET_COMPLEX (&z, w, ai / (2.0 * w));
  302.     }
  303.       else
  304.     {
  305.       double ai = GSL_IMAG (a);
  306.       double vi = (ai >= 0) ? w : -w;
  307.       GSL_SET_COMPLEX (&z, ai / (2.0 * vi), vi);
  308.     }
  309.     }
  310.  
  311.   return z;
  312. }
  313.  
  314. gsl_complex
  315. gsl_complex_sqrt_real (double x)
  316. {                /* z=sqrt(x) */
  317.   gsl_complex z;
  318.  
  319.   if (x >= 0)
  320.     {
  321.       GSL_SET_COMPLEX (&z, sqrt (x), 0.0);
  322.     }
  323.   else
  324.     {
  325.       GSL_SET_COMPLEX (&z, 0.0, sqrt (-x));
  326.     }
  327.  
  328.   return z;
  329. }
  330.  
  331. gsl_complex
  332. gsl_complex_exp (gsl_complex a)
  333. {                /* z=exp(a) */
  334.   double rho = exp (GSL_REAL (a));
  335.   double theta = GSL_IMAG (a);
  336.  
  337.   gsl_complex z;
  338.   GSL_SET_COMPLEX (&z, rho * cos (theta), rho * sin (theta));
  339.   return z;
  340. }
  341.  
  342. gsl_complex
  343. gsl_complex_pow (gsl_complex a, gsl_complex b)
  344. {                /* z=a^b */
  345.   gsl_complex z;
  346.  
  347.   if (GSL_REAL (a) == 0 && GSL_IMAG (a) == 0.0)
  348.     {
  349.       GSL_SET_COMPLEX (&z, 0.0, 0.0);
  350.     }
  351.   else
  352.     {
  353.       double logr = gsl_complex_logabs (a);
  354.       double theta = gsl_complex_arg (a);
  355.  
  356.       double br = GSL_REAL (b), bi = GSL_IMAG (b);
  357.  
  358.       double rho = exp (logr * br - bi * theta);
  359.       double beta = theta * br + bi * logr;
  360.  
  361.       GSL_SET_COMPLEX (&z, rho * cos (beta), rho * sin (beta));
  362.     }
  363.  
  364.   return z;
  365. }
  366.  
  367. gsl_complex
  368. gsl_complex_pow_real (gsl_complex a, double b)
  369. {                /* z=a^b */
  370.   gsl_complex z;
  371.  
  372.   if (GSL_REAL (a) == 0 && GSL_IMAG (a) == 0)
  373.     {
  374.       GSL_SET_COMPLEX (&z, 0, 0);
  375.     }
  376.   else
  377.     {
  378.       double logr = gsl_complex_logabs (a);
  379.       double theta = gsl_complex_arg (a);
  380.       double rho = exp (logr * b);
  381.       double beta = theta * b;
  382.       GSL_SET_COMPLEX (&z, rho * cos (beta), rho * sin (beta));
  383.     }
  384.  
  385.   return z;
  386. }
  387.  
  388. gsl_complex
  389. gsl_complex_log (gsl_complex a)
  390. {                /* z=log(a) */
  391.   double logr = gsl_complex_logabs (a);
  392.   double theta = gsl_complex_arg (a);
  393.  
  394.   gsl_complex z;
  395.   GSL_SET_COMPLEX (&z, logr, theta);
  396.   return z;
  397. }
  398.  
  399. gsl_complex
  400. gsl_complex_log10 (gsl_complex a)
  401. {                /* z = log10(a) */
  402.   return gsl_complex_mul_real (gsl_complex_log (a), 1 / log (10.));
  403. }
  404.  
  405. gsl_complex
  406. gsl_complex_log_b (gsl_complex a, gsl_complex b)
  407. {
  408.   return gsl_complex_div (gsl_complex_log (a), gsl_complex_log (b));
  409. }
  410.  
  411. /***********************************************************************
  412.  * Complex trigonometric functions
  413.  ***********************************************************************/
  414.  
  415. gsl_complex
  416. gsl_complex_sin (gsl_complex a)
  417. {                /* z = sin(a) */
  418.   double R = GSL_REAL (a), I = GSL_IMAG (a);
  419.  
  420.   gsl_complex z;
  421.  
  422.   if (I == 0.0) 
  423.     {
  424.       /* avoid returing negative zero (-0.0) for the imaginary part  */
  425.  
  426.       GSL_SET_COMPLEX (&z, sin (R), 0.0);  
  427.     } 
  428.   else 
  429.     {
  430.       GSL_SET_COMPLEX (&z, sin (R) * cosh (I), cos (R) * sinh (I));
  431.     }
  432.  
  433.   return z;
  434. }
  435.  
  436. gsl_complex
  437. gsl_complex_cos (gsl_complex a)
  438. {                /* z = cos(a) */
  439.   double R = GSL_REAL (a), I = GSL_IMAG (a);
  440.  
  441.   gsl_complex z;
  442.  
  443.   if (I == 0.0) 
  444.     {
  445.       /* avoid returing negative zero (-0.0) for the imaginary part  */
  446.  
  447.       GSL_SET_COMPLEX (&z, cos (R), 0.0);  
  448.     } 
  449.   else 
  450.     {
  451.       GSL_SET_COMPLEX (&z, cos (R) * cosh (I), sin (R) * sinh (-I));
  452.     }
  453.  
  454.   return z;
  455. }
  456.  
  457. gsl_complex
  458. gsl_complex_tan (gsl_complex a)
  459. {                /* z = tan(a) */
  460.   double R = GSL_REAL (a), I = GSL_IMAG (a);
  461.  
  462.   gsl_complex z;
  463.  
  464.   if (fabs (I) < 1)
  465.     {
  466.       double D = pow (cos (R), 2.0) + pow (sinh (I), 2.0);
  467.  
  468.       GSL_SET_COMPLEX (&z, 0.5 * sin (2 * R) / D, 0.5 * sinh (2 * I) / D);
  469.     }
  470.   else
  471.     {
  472.       double u = exp (-I);
  473.       double C = 2 * u / (1 - pow (u, 2.0));
  474.       double D = 1 + pow (cos (R), 2.0) * pow (C, 2.0);
  475.  
  476.       double S = pow (C, 2.0);
  477.       double T = 1.0 / tanh (I);
  478.  
  479.       GSL_SET_COMPLEX (&z, 0.5 * sin (2 * R) * S / D, T / D);
  480.     }
  481.  
  482.   return z;
  483. }
  484.  
  485. gsl_complex
  486. gsl_complex_sec (gsl_complex a)
  487. {                /* z = sec(a) */
  488.   gsl_complex z = gsl_complex_cos (a);
  489.   return gsl_complex_inverse (z);
  490. }
  491.  
  492. gsl_complex
  493. gsl_complex_csc (gsl_complex a)
  494. {                /* z = csc(a) */
  495.   gsl_complex z = gsl_complex_sin (a);
  496.   return gsl_complex_inverse(z);
  497. }
  498.  
  499.  
  500. gsl_complex
  501. gsl_complex_cot (gsl_complex a)
  502. {                /* z = cot(a) */
  503.   gsl_complex z = gsl_complex_tan (a);
  504.   return gsl_complex_inverse (z);
  505. }
  506.  
  507. /**********************************************************************
  508.  * Inverse Complex Trigonometric Functions
  509.  **********************************************************************/
  510.  
  511. gsl_complex
  512. gsl_complex_arcsin (gsl_complex a)
  513. {                /* z = arcsin(a) */
  514.   double R = GSL_REAL (a), I = GSL_IMAG (a);
  515.   gsl_complex z;
  516.  
  517.   if (I == 0)
  518.     {
  519.       z = gsl_complex_arcsin_real (R);
  520.     }
  521.   else
  522.     {
  523.       double x = fabs (R), y = fabs (I);
  524.       double r = hypot (x + 1, y), s = hypot (x - 1, y);
  525.       double A = 0.5 * (r + s);
  526.       double B = x / A;
  527.       double y2 = y * y;
  528.  
  529.       double real, imag;
  530.  
  531.       const double A_crossover = 1.5, B_crossover = 0.6417;
  532.  
  533.       if (B <= B_crossover)
  534.     {
  535.       real = asin (B);
  536.     }
  537.       else
  538.     {
  539.       if (x <= 1)
  540.         {
  541.           double D = 0.5 * (A + x) * (y2 / (r + x + 1) + (s + (1 - x)));
  542.           real = atan (x / sqrt (D));
  543.         }
  544.       else
  545.         {
  546.           double Apx = A + x;
  547.           double D = 0.5 * (Apx / (r + x + 1) + Apx / (s + (x - 1)));
  548.           real = atan (x / (y * sqrt (D)));
  549.         }
  550.     }
  551.  
  552.       if (A <= A_crossover)
  553.     {
  554.       double Am1;
  555.  
  556.       if (x < 1)
  557.         {
  558.           Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 / (s + (1 - x)));
  559.         }
  560.       else
  561.         {
  562.           Am1 = 0.5 * (y2 / (r + (x + 1)) + (s + (x - 1)));
  563.         }
  564.  
  565.       imag = log1p (Am1 + sqrt (Am1 * (A + 1)));
  566.     }
  567.       else
  568.     {
  569.       imag = log (A + sqrt (A * A - 1));
  570.     }
  571.  
  572.       GSL_SET_COMPLEX (&z, (R >= 0) ? real : -real, (I >= 0) ? imag : -imag);
  573.     }
  574.  
  575.   return z;
  576. }
  577.  
  578. gsl_complex
  579. gsl_complex_arcsin_real (double a)
  580. {                /* z = arcsin(a) */
  581.   gsl_complex z;
  582.  
  583.   if (fabs (a) <= 1.0)
  584.     {
  585.       GSL_SET_COMPLEX (&z, asin (a), 0.0);
  586.     }
  587.   else
  588.     {
  589.       if (a < 0.0)
  590.     {
  591.       GSL_SET_COMPLEX (&z, -M_PI_2, acosh (-a));
  592.     }
  593.       else
  594.     {
  595.       GSL_SET_COMPLEX (&z, M_PI_2, -acosh (a));
  596.     }
  597.     }
  598.  
  599.   return z;
  600. }
  601.  
  602. gsl_complex
  603. gsl_complex_arccos (gsl_complex a)
  604. {                /* z = arccos(a) */
  605.   double R = GSL_REAL (a), I = GSL_IMAG (a);
  606.   gsl_complex z;
  607.  
  608.   if (I == 0)
  609.     {
  610.       z = gsl_complex_arccos_real (R);
  611.     }
  612.   else
  613.     {
  614.       double x = fabs (R), y = fabs (I);
  615.       double r = hypot (x + 1, y), s = hypot (x - 1, y);
  616.       double A = 0.5 * (r + s);
  617.       double B = x / A;
  618.       double y2 = y * y;
  619.  
  620.       double real, imag;
  621.  
  622.       const double A_crossover = 1.5, B_crossover = 0.6417;
  623.  
  624.       if (B <= B_crossover)
  625.     {
  626.       real = acos (B);
  627.     }
  628.       else
  629.     {
  630.       if (x <= 1)
  631.         {
  632.           double D = 0.5 * (A + x) * (y2 / (r + x + 1) + (s + (1 - x)));
  633.           real = atan (sqrt (D) / x);
  634.         }
  635.       else
  636.         {
  637.           double Apx = A + x;
  638.           double D = 0.5 * (Apx / (r + x + 1) + Apx / (s + (x - 1)));
  639.           real = atan ((y * sqrt (D)) / x);
  640.         }
  641.     }
  642.  
  643.       if (A <= A_crossover)
  644.     {
  645.       double Am1;
  646.  
  647.       if (x < 1)
  648.         {
  649.           Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 / (s + (1 - x)));
  650.         }
  651.       else
  652.         {
  653.           Am1 = 0.5 * (y2 / (r + (x + 1)) + (s + (x - 1)));
  654.         }
  655.  
  656.       imag = log1p (Am1 + sqrt (Am1 * (A + 1)));
  657.     }
  658.       else
  659.     {
  660.       imag = log (A + sqrt (A * A - 1));
  661.     }
  662.  
  663.       GSL_SET_COMPLEX (&z, (R >= 0) ? real : M_PI - real, (I >= 0) ? -imag : imag);
  664.     }
  665.  
  666.   return z;
  667. }
  668.  
  669. gsl_complex
  670. gsl_complex_arccos_real (double a)
  671. {                /* z = arccos(a) */
  672.   gsl_complex z;
  673.  
  674.   if (fabs (a) <= 1.0)
  675.     {
  676.       GSL_SET_COMPLEX (&z, acos (a), 0);
  677.     }
  678.   else
  679.     {
  680.       if (a < 0.0)
  681.     {
  682.       GSL_SET_COMPLEX (&z, M_PI, -acosh (-a));
  683.     }
  684.       else
  685.     {
  686.       GSL_SET_COMPLEX (&z, 0, acosh (a));
  687.     }
  688.     }
  689.  
  690.   return z;
  691. }
  692.  
  693. gsl_complex
  694. gsl_complex_arctan (gsl_complex a)
  695. {                /* z = arctan(a) */
  696.   double R = GSL_REAL (a), I = GSL_IMAG (a);
  697.   gsl_complex z;
  698.  
  699.   if (I == 0)
  700.     {
  701.       GSL_SET_COMPLEX (&z, atan (R), 0);
  702.     }
  703.   else
  704.     {
  705.       /* FIXME: This is a naive implementation which does not fully
  706.          take into account cancellation errors, overflow, underflow
  707.          etc.  It would benefit from the Hull et al treatment. */
  708.  
  709.       double r = hypot (R, I);
  710.  
  711.       double imag;
  712.  
  713.       double u = 2 * I / (1 + r * r);
  714.  
  715.       /* FIXME: the following cross-over should be optimized but 0.1
  716.          seems to work ok */
  717.  
  718.       if (fabs (u) < 0.1)
  719.     {
  720.       imag = 0.25 * (log1p (u) - log1p (-u));
  721.     }
  722.       else
  723.     {
  724.       double A = hypot (R, I + 1);
  725.       double B = hypot (R, I - 1);
  726.       imag = 0.5 * log (A / B);
  727.     }
  728.  
  729.       if (R == 0)
  730.     {
  731.       if (I > 1)
  732.         {
  733.           GSL_SET_COMPLEX (&z, M_PI_2, imag);
  734.         }
  735.       else if (I < -1)
  736.         {
  737.           GSL_SET_COMPLEX (&z, -M_PI_2, imag);
  738.         }
  739.       else
  740.         {
  741.           GSL_SET_COMPLEX (&z, 0, imag);
  742.         };
  743.     }
  744.       else
  745.     {
  746.       GSL_SET_COMPLEX (&z, 0.5 * atan2 (2 * R, ((1 + r) * (1 - r))), imag);
  747.     }
  748.     }
  749.  
  750.   return z;
  751. }
  752.  
  753. gsl_complex
  754. gsl_complex_arcsec (gsl_complex a)
  755. {                /* z = arcsec(a) */
  756.   gsl_complex z = gsl_complex_inverse (a);
  757.   return gsl_complex_arccos (z);
  758. }
  759.  
  760. gsl_complex
  761. gsl_complex_arcsec_real (double a)
  762. {                /* z = arcsec(a) */
  763.   gsl_complex z;
  764.  
  765.   if (a <= -1.0 || a >= 1.0)
  766.     {
  767.       GSL_SET_COMPLEX (&z, acos (1 / a), 0.0);
  768.     }
  769.   else
  770.     {
  771.       if (a >= 0.0)
  772.     {
  773.       GSL_SET_COMPLEX (&z, 0, acosh (1 / a));
  774.     }
  775.       else
  776.     {
  777.       GSL_SET_COMPLEX (&z, M_PI, -acosh (-1 / a));
  778.     }
  779.     }
  780.  
  781.   return z;
  782. }
  783.  
  784. gsl_complex
  785. gsl_complex_arccsc (gsl_complex a)
  786. {                /* z = arccsc(a) */
  787.   gsl_complex z = gsl_complex_inverse (a);
  788.   return gsl_complex_arcsin (z);
  789. }
  790.  
  791. gsl_complex
  792. gsl_complex_arccsc_real (double a)
  793. {                /* z = arccsc(a) */
  794.   gsl_complex z;
  795.  
  796.   if (a <= -1.0 || a >= 1.0)
  797.     {
  798.       GSL_SET_COMPLEX (&z, asin (1 / a), 0.0);
  799.     }
  800.   else
  801.     {
  802.       if (a >= 0.0)
  803.     {
  804.       GSL_SET_COMPLEX (&z, M_PI_2, -acosh (1 / a));
  805.     }
  806.       else
  807.     {
  808.       GSL_SET_COMPLEX (&z, -M_PI_2, -acosh (-1 / a));
  809.     }
  810.     }
  811.  
  812.   return z;
  813. }
  814.  
  815. gsl_complex
  816. gsl_complex_arccot (gsl_complex a)
  817. {                /* z = arccot(a) */
  818.   gsl_complex z;
  819.  
  820.   if (GSL_REAL (a) == 0.0 && GSL_IMAG (a) == 0.0)
  821.     {
  822.       GSL_SET_COMPLEX (&z, M_PI_2, 0);
  823.     }
  824.   else
  825.     {
  826.       z = gsl_complex_inverse (a);
  827.       z = gsl_complex_arctan (z);
  828.     }
  829.  
  830.   return z;
  831. }
  832.  
  833. /**********************************************************************
  834.  * Complex Hyperbolic Functions
  835.  **********************************************************************/
  836.  
  837. gsl_complex
  838. gsl_complex_sinh (gsl_complex a)
  839. {                /* z = sinh(a) */
  840.   double R = GSL_REAL (a), I = GSL_IMAG (a);
  841.  
  842.   gsl_complex z;
  843.   GSL_SET_COMPLEX (&z, sinh (R) * cos (I), cosh (R) * sin (I));
  844.   return z;
  845. }
  846.  
  847. gsl_complex
  848. gsl_complex_cosh (gsl_complex a)
  849. {                /* z = cosh(a) */
  850.   double R = GSL_REAL (a), I = GSL_IMAG (a);
  851.  
  852.   gsl_complex z;
  853.   GSL_SET_COMPLEX (&z, cosh (R) * cos (I), sinh (R) * sin (I));
  854.   return z;
  855. }
  856.  
  857. gsl_complex
  858. gsl_complex_tanh (gsl_complex a)
  859. {                /* z = tanh(a) */
  860.   double R = GSL_REAL (a), I = GSL_IMAG (a);
  861.  
  862.   gsl_complex z;
  863.  
  864.   if (fabs(R) < 1.0) 
  865.     {
  866.       double D = pow (cos (I), 2.0) + pow (sinh (R), 2.0);
  867.       
  868.       GSL_SET_COMPLEX (&z, sinh (R) * cosh (R) / D, 0.5 * sin (2 * I) / D);
  869.     }
  870.   else
  871.     {
  872.       double D = pow (cos (I), 2.0) + pow (sinh (R), 2.0);
  873.       double F = 1 + pow (cos (I) / sinh (R), 2.0);
  874.  
  875.       GSL_SET_COMPLEX (&z, 1.0 / (tanh (R) * F), 0.5 * sin (2 * I) / D);
  876.     }
  877.  
  878.   return z;
  879. }
  880.  
  881. gsl_complex
  882. gsl_complex_sech (gsl_complex a)
  883. {                /* z = sech(a) */
  884.   gsl_complex z = gsl_complex_cosh (a);
  885.   return gsl_complex_inverse (z);
  886. }
  887.  
  888. gsl_complex
  889. gsl_complex_csch (gsl_complex a)
  890. {                /* z = csch(a) */
  891.   gsl_complex z = gsl_complex_sinh (a);
  892.   return gsl_complex_inverse (z);
  893. }
  894.  
  895. gsl_complex
  896. gsl_complex_coth (gsl_complex a)
  897. {                /* z = coth(a) */
  898.   gsl_complex z = gsl_complex_tanh (a);
  899.   return gsl_complex_inverse (z);
  900. }
  901.  
  902. /**********************************************************************
  903.  * Inverse Complex Hyperbolic Functions
  904.  **********************************************************************/
  905.  
  906. gsl_complex
  907. gsl_complex_arcsinh (gsl_complex a)
  908. {                /* z = arcsinh(a) */
  909.   gsl_complex z = gsl_complex_mul_imag(a, 1.0);
  910.   z = gsl_complex_arcsin (z);
  911.   z = gsl_complex_mul_imag (z, -1.0);
  912.   return z;
  913. }
  914.  
  915. gsl_complex
  916. gsl_complex_arccosh (gsl_complex a)
  917. {                /* z = arccosh(a) */
  918.   gsl_complex z = gsl_complex_arccos (a);
  919.   z = gsl_complex_mul_imag (z, GSL_IMAG(z) > 0 ? -1.0 : 1.0);
  920.   return z;
  921. }
  922.  
  923. gsl_complex
  924. gsl_complex_arccosh_real (double a)
  925. {                /* z = arccosh(a) */
  926.   gsl_complex z;
  927.  
  928.   if (a >= 1)
  929.     {
  930.       GSL_SET_COMPLEX (&z, acosh (a), 0);
  931.     }
  932.   else
  933.     {
  934.       if (a >= -1.0)
  935.     {
  936.       GSL_SET_COMPLEX (&z, 0, acos (a));
  937.     }
  938.       else
  939.     {
  940.       GSL_SET_COMPLEX (&z, acosh (-a), M_PI);
  941.     }
  942.     }
  943.  
  944.   return z;
  945. }
  946.  
  947. gsl_complex
  948. gsl_complex_arctanh (gsl_complex a)
  949. {                /* z = arctanh(a) */
  950.   if (GSL_IMAG (a) == 0.0)
  951.     {
  952.       return gsl_complex_arctanh_real (GSL_REAL (a));
  953.     }
  954.   else
  955.     {
  956.       gsl_complex z = gsl_complex_mul_imag(a, 1.0);
  957.       z = gsl_complex_arctan (z);
  958.       z = gsl_complex_mul_imag (z, -1.0);
  959.       return z;
  960.     }
  961. }
  962.  
  963. gsl_complex
  964. gsl_complex_arctanh_real (double a)
  965. {                /* z = arctanh(a) */
  966.   gsl_complex z;
  967.  
  968.   if (a > -1.0 && a < 1.0)
  969.     {
  970.       GSL_SET_COMPLEX (&z, atanh (a), 0);
  971.     }
  972.   else
  973.     {
  974.       GSL_SET_COMPLEX (&z, atanh (1 / a), (a < 0) ? M_PI_2 : -M_PI_2);
  975.     }
  976.  
  977.   return z;
  978. }
  979.  
  980. gsl_complex
  981. gsl_complex_arcsech (gsl_complex a)
  982. {                /* z = arcsech(a); */
  983.   gsl_complex t = gsl_complex_inverse (a);
  984.   return gsl_complex_arccosh (t);
  985. }
  986.  
  987. gsl_complex
  988. gsl_complex_arccsch (gsl_complex a)
  989. {                /* z = arccsch(a) */
  990.   gsl_complex t = gsl_complex_inverse (a);
  991.   return gsl_complex_arcsinh (t);
  992. }
  993.  
  994. gsl_complex
  995. gsl_complex_arccoth (gsl_complex a)
  996. {                /* z = arccoth(a) */
  997.   gsl_complex t = gsl_complex_inverse (a);
  998.   return gsl_complex_arctanh (t);
  999. }
  1000.